%% ARIMA模型预测ARIMA model prediction
clear,clc
%% 原始数据列：
n=1000;
x=cumsum(ones(1,n));
y=1/n*x+sin(x/(4*pi))/10+randn(1,n);% 具有上升趋势，小幅周期波动,噪声的数据列
%% 判断数据列平稳性并做差分处理
dm=4;% 约束最大进行4次差分
yd=y;% 差分数列初始化
for dp=1:dm+2
    n1=length(yd);
    ry=sort(yd);Ry=zeros(1,n1);
    for i=1:n1
        [~,r]=find(yd==ry(i));
        Ry(i)=r;% 获得y数据列秩统计量
    end
    xp=cumsum(ones(1,n1));
    d=xp-Ry;
    Qxy=1-6/(n1*(n1^2-1))*sum(d.^2);
    T=Qxy*sqrt(n1-2)/sqrt(1-Qxy^2);
    t_0=tinv(0.975,n1-2);% 计算上alpha/2分位数
    if(abs(T)<=t_0)
        break
    elseif(dp==dm+2)
        disp('原始数据列无法差分平稳化，可能不适合ARIMA模型')
        break
    else
        yd=diff(yd);
    end
end
d_p=dp-1;% 最终差分次数
%% 自相关系数，偏自相关系数计算
y_mean=mean(yd);% 数据yd平均值
y_0=yd-y_mean;% 序列减去其均值
rho=zeros(1,n1);
for k=0:n1-1
    rho(k+1)=sum(y_0(1:n1-k).*y_0(k+1:n1))/sum(y_0.^2);
end
phi=zeros(1,n1-1);
for k=0:n1-2
    tic
    rho1=rho(2:k+2);
    kp=length(rho1);
    rho2=zeros(kp);
    for k0=1:kp
        rho2(k0,k0:end)=rho(1:kp-k0+1);
        if(k0>=2)
            rho2(k0,1:k0-1)=fliplr(rho1(1:k0-1));
        end
    end
    phi1=rho2\rho1';phi(k+1)=phi1(end);
    toc
end
%% 定阶预测
% ARMA(p,q)预测一般需满足：数据自相关系数(ACF)q阶后衰减趋于0，偏自相关系数
% (PACF)p阶后衰减趋于0.
% 如果样本ACF和样本PACF在最初k阶明显大于2倍标准差，而后几乎95%的系数都落在2
% 倍标准差的范围内，且非零系数衰减为小值波动的过程非常突然，通常视为k阶截尾；
% 如果有超过5%的样本相关系数大于2倍标准差，或者非零系数衰减为小值波动的过程
% 比较缓慢或连续，通常视为拖尾。
rho_std=std(rho);% 自相关系数标准差
phi_std=std(phi);% 偏自相关系数标准差
p=0;q=0;
for i=rho
    if(abs(i)>2*rho_std)
        p=p+1;
    else
        break
    end
end
for i=phi
    if(abs(i)>2*phi_std)
        q=q+1;
    else
        break
    end
end
% 以上程序简单进行了定阶，通常利用贝叶斯信息等可以更准确的进行定阶，限于篇幅
% 本文不再赘述。
% 利用matlab内置函数进行预测；
AR_Order=p;
MA_Order=q;
AR=arima(AR_Order, d_p, MA_Order);
EstMdl=estimate(AR,y');
step=100;% 预测步数
fore_y=forecast(EstMdl,step,'Y0',y');
figure(1)
plot([y,fore_y'],'r*-')
hold on
plot(y,'b*-')
hold off
legend('预测值','实际值')